version 16
cap mata mata drop _ddf()
cap mata mata drop _nddf()
cap mata mata drop _ddf2()
cap mata mata drop _nddf2()
cap mata mata drop _tenddf()
cap mata mata drop _tenddf2()
mata:
mata clear 








void function _ddf( string scalar     d, ///
                    real scalar       k1, ///
                    real scalar       k2, ///
					string scalar     touse1,  ///
					string scalar     touse2, ///
					string scalar     gname, ///
					string scalar     bname, ///
					real scalar       rstype, ///					
					real scalar       maxiter, ///
					real scalar       tol)
	{


          data=st_data(.,d,touse1)
          data=data'
          M=rows(data)          
          g=st_data(.,gname,touse1)
          g=g'
        
          dataref=st_data(.,d,touse2)
          dataref=dataref'

          Xref=dataref[1..k1,.]
          Yref=dataref[k1+1..(k1+k2),.]
          Bref=dataref[(k1+k2+1)..M,.]
          X=data[1..k1,.]
          Y=data[k1+1..(k1+k2),.]
          B=data[(k1+k2+1)..M,.]
          theta=J(cols(data),1,.)         
		 
		 for(j=1;j<=cols(X);j++){
			theta[j]=_ddf2(X[.,j],Y[.,j],B[.,j],Xref,Yref,Bref,g[.,j],rstype,maxiter,tol)
		 }

          st_view(gen=.,.,bname,touse1)
          gen[.,.]=theta
		

}




									
									
real matrix function ddfcontemporaneous(    real matrix       data0, ///
											real matrix       dataref, ///
											real scalar       k1, ///
											real scalar       k2, ///
											real matrix       gmat, ///
											real scalar       rstype, ///					
											real scalar       maxiter, ///
											real scalar       tol)									
	{


		  t0  = data0[.,2]
		  t1  = dataref[.,2]
          Mc=cols(data) 
		  
          
	      t0uniq=uniqrows(t0) 
		  D0=J(rows(data0),1,.)  
		  for(j=1;j<=length(t0uniq);j++){
		      XYB    = select(data0[.,2...Mc],t0:==t0uniq[j])
			  XYBref = select(dataref[.,2...Mc],t1:==t0uniq[j])
			  gmat2  = select(gmat,t0:==t0uniq[j])
			  pos=select(1::rows(data0),t0:==t0uniq[j])
			  D0[pos,1]= calddf(XYB',XYBref',k1,k2,gmat2',rstype,maxiter,tol)
		  }
		  
		  return(D0)


}



real matrix function ddfcontemporaneousbt(  real matrix       data0, ///
											real matrix       dataref, ///
											real scalar       k1, ///
											real scalar       k2, ///
											real matrix       gmat, ///
										    real scalar       gamma, ///
										    real scalar       level, ///											
											real scalar       brep,
											real scalar       rstype, ///					
											real scalar       maxiter, ///
											real scalar       tol)									
	{
         // sort by id t
		 
		  struct bootdea scalar cires
          msub=ceil(length(uniqrows(dataref[.,1]))^gamma)
		  t0  = data0[.,2]
		  t1  = dataref[.,2]
          Mc=cols(data) 
	      t0uniq=uniqrows(t0) 
		  //Kr=length(uniqrows(dataref[.,1]))
		  D0=J(rows(data0),1,.)  
		  for(j=1;j<=length(t0uniq);j++){
		      XYB    = select(data0[.,2...Mc],t0:==t0uniq[j])
			  XYBref = select(dataref[.,2...Mc],t1:==t0uniq[j])
			  gmat2  = select(gmat,t0:==t0uniq[j])
			  pos=select(1::rows(data0),t0:==t0uniq[j])
			  D0[pos,1]= calddf(XYB',XYBref',k1,k2,gmat2',rstype,maxiter,tol)
		  }
		  
		  id = dataref[.,1]
		  mm_panels(id,idinfo=.)
		  
		  Dboot=J(rows(data0),brep,.)
		  
		  for(i=1;i<=brep;i++){
		      idx=mm_sample(msub,.,idinfo,.,1)
			  datastar=dataref[idx,.]
			  Dboot[.,i]=ddfcontemporaneous(data0,datastar,k1,k2,gmat,rstype,maxiter,tol)
		  }
		  
		  
		 for(j=1;j<=length(t0uniq);j++){
			  pos=select(1::rows(data0),t0:==t0uniq[j])
			  Kr=length(select(dataref[.,1],t1:==t0uniq[j]))
			  cires=biasandciddf(D0[pos,1],Dboot[pos,.], k1, cols(data)-k1-2, length(pos), Kr, level,0,msub)
			  Dbc[pos,1]=cires.tebc
			  LL[pos,1] =cires.LL
			  UU[pos,1] =cires.UU
			  bias[pos,1]=cires.bias
			  vari[pos,1]=cires.vari
		  }
		  


}




real matrix function ddfseq(                real matrix       data0, ///
											real matrix       dataref, ///
											real scalar       k1, ///
											real scalar       k2, ///
											real matrix       gmat, ///
											real scalar       rstype, ///					
											real scalar       maxiter, ///
											real scalar       tol)									
	{


		  t0  = data0[.,2]
		  t1  = dataref[.,2]
          Mc=cols(data) 
		  
          
	      t0uniq=uniqrows(t0) 
		  D0=J(rows(data0),1,.)  
		  for(j=1;j<=length(t0uniq);j++){
		      XYB    = select(data0[.,2...Mc],t0:==t0uniq[j])
			  XYBref = select(dataref[.,2...Mc],t1:<=t0uniq[j])
			  gmat2  = select(gmat,t0:==t0uniq[j])
			  pos=select(1::rows(data0),t0:==t0uniq[j])
			  D0[pos,1]= calddf(XYB',XYBref',k1,k2,gmat2',rstype,maxiter,tol)
		  }
		  
		  return(D0)


}



real matrix function ddfseqbt(              real matrix       data0, ///
											real matrix       dataref, ///
											real scalar       k1, ///
											real scalar       k2, ///
											real matrix       gmat, ///
										    real scalar       gamma, ///
										    real scalar       level, ///											
											real scalar       brep,
											real scalar       rstype, ///					
											real scalar       maxiter, ///
											real scalar       tol)									
	{
         // sort by id t
		 
		  struct bootdea scalar cires
          msub=ceil(length(uniqrows(dataref[.,1]))^gamma)
		  t0  = data0[.,2]
		  t1  = dataref[.,2]
          Mc=cols(data) 
	      t0uniq=uniqrows(t0) 
		  //Kr=length(uniqrows(dataref[.,1]))
		  D0=J(rows(data0),1,.)  
		  for(j=1;j<=length(t0uniq);j++){
		      XYB    = select(data0[.,2...Mc],t0:==t0uniq[j])
			  XYBref = select(dataref[.,2...Mc],t1:<=t0uniq[j])
			  gmat2  = select(gmat,t0:==t0uniq[j])
			  pos=select(1::rows(data0),t0:==t0uniq[j])
			  D0[pos,1]= calddf(XYB',XYBref',k1,k2,gmat2',rstype,maxiter,tol)
		  }
		  
		  id = dataref[.,1]
		  mm_panels(id,idinfo=.)
		  
		  Dboot=J(rows(data0),brep,.)
		  
		  for(i=1;i<=brep;i++){
		      idx=mm_sample(msub,.,idinfo,.,1)
			  datastar=dataref[idx,.]
			  Dboot[.,i]=ddfseq(data0,datastar,k1,k2,gmat,rstype,maxiter,tol)
		  }
		  
		  
		 for(j=1;j<=length(t0uniq);j++){
			  pos=select(1::rows(data0),t0:==t0uniq[j])
			  Kr=length(select(dataref[.,1],t1:<=t0uniq[j]))
			  cires=biasandciddf(D0[pos,1],Dboot[pos,.], k1, cols(data)-k1-2, length(pos), Kr, level,0,msub*sum(t1:<=t0uniq[j]))
			  Dbc[pos,1]=cires.tebc
			  LL[pos,1] =cires.LL
			  UU[pos,1] =cires.UU
			  bias[pos,1]=cires.bias
			  vari[pos,1]=cires.vari
		  }
		  


}




real matrix function ddfglobal(                real matrix       data0, ///
											real matrix       dataref, ///
											real scalar       k1, ///
											real scalar       k2, ///
											real matrix       gmat, ///
											real scalar       rstype, ///					
											real scalar       maxiter, ///
											real scalar       tol)									
{


		XYB    = data0[.,2...cols(data0)]
	    XYBref = dataref[.,2...cols(dataref)]
		D0 = calddf(XYB',XYBref',k1,k2,gmat',rstype,maxiter,tol)
		return(D0)


}



real matrix function ddfglobalbt(              real matrix       data0, ///
											real matrix       dataref, ///
											real scalar       k1, ///
											real scalar       k2, ///
											real matrix       gmat, ///
										    real scalar       gamma, ///
										    real scalar       level, ///											
											real scalar       brep,
											real scalar       rstype, ///					
											real scalar       maxiter, ///
											real scalar       tol)									
	{
         // sort by id t
		 
		  struct bootdea scalar cires
          msub=ceil(length(uniqrows(dataref[.,1]))^gamma)
		  t0  = data0[.,2]
		  t1  = dataref[.,2]
          Mc=cols(data) 
	      t0uniq=uniqrows(t0) 
		  Kr=rows(dataref)
		
		  
		XYB    = data0[.,2...cols(data0)]
		XYBref = dataref[.,2...cols(dataref)]
		D0 = calddf(XYB',XYBref',k1,k2,gmat',rstype,maxiter,tol)
		  
		  
		  id = dataref[.,1]
		  mm_panels(id,idinfo=.)
		  
		  Dboot=J(rows(data0),brep,.)
		  msubmat=J(1,brep,.)
		  for(i=1;i<=brep;i++){
		      idx=mm_sample(msub,.,idinfo,.,1)
			  datastar=dataref[idx,.]
			  msubmat[1,i]=rows(datastar)
			  Dboot[.,i]=ddfglobal(data0,datastar,k1,k2,gmat,rstype,maxiter,tol)
		  }
		  
		  
		 for(j=1;j<=length(t0uniq);j++){
			  pos=select(1::rows(data0),t0:==t0uniq[j])
			  cires=biasandciddf(D0[pos,1],Dboot[pos,.], k1, cols(data)-k1-2, length(pos), Kr, level,0,msubmat[1,j])
			  Dbc[pos,1]=cires.tebc
			  LL[pos,1] =cires.LL
			  UU[pos,1] =cires.UU
			  bias[pos,1]=cires.bias
			  vari[pos,1]=cires.vari
		  }
		  


}







real matrix function ddfwindow(             real matrix       data0, ///
											real matrix       dataref, ///
											real scalar       band, ///
											real scalar       k1, ///
											real scalar       k2, ///
											real matrix       gmat, ///
											real scalar       rstype, ///					
											real scalar       maxiter, ///
											real scalar       tol)									
	{


		  t0  = data0[.,2]
		  t1  = dataref[.,2]
          Mc=cols(data) 
		  
          
	      t0uniq=uniqrows(t0) 
		  D0=J(rows(data0),1,.)  
		  for(j=1;j<=length(t0uniq);j++){
		      XYB    = select(data0[.,2...Mc],t0:==t0uniq[j])
			  XYBref = select(dataref[.,2...Mc],t1:<=t0uniq[j]+band & t1:>=t0uniq[j]-band)
			  gmat2  = select(gmat,t0:==t0uniq[j])
			  pos=select(1::rows(data0),t0:==t0uniq[j])
			  D0[pos,1]= calddf(XYB',XYBref',k1,k2,gmat2',rstype,maxiter,tol)
		  }
		  
		  return(D0)


}



real matrix function ddfwindowbt(              real matrix       data0, ///
											real matrix       dataref, ///
											real scalar       band,
											real scalar       k1, ///
											real scalar       k2, ///
											real matrix       gmat, ///
										    real scalar       gamma, ///
										    real scalar       level, ///											
											real scalar       brep,
											real scalar       rstype, ///					
											real scalar       maxiter, ///
											real scalar       tol)									
	{
         // sort by id t
		 
		  struct bootdea scalar cires
          msub=ceil(length(uniqrows(dataref[.,1]))^gamma)
		  t0  = data0[.,2]
		  t1  = dataref[.,2]
          Mc=cols(data) 
	      t0uniq=uniqrows(t0) 
		  //Kr=length(uniqrows(dataref[.,1]))
		  D0=J(rows(data0),1,.)  
		  for(j=1;j<=length(t0uniq);j++){
		      XYB    = select(data0[.,2...Mc],t0:==t0uniq[j])
			  XYBref = select(dataref[.,2...Mc],t1:<=t0uniq[j]+band & t1:>=t0uniq[j]-band)
			  gmat2  = select(gmat,t0:==t0uniq[j])
			  pos=select(1::rows(data0),t0:==t0uniq[j])
			  D0[pos,1]= calddf(XYB',XYBref',k1,k2,gmat2',rstype,maxiter,tol)
		  }
		  
		  id = dataref[.,1]
		  mm_panels(id,idinfo=.)
		  
		  Dboot=J(rows(data0),brep,.)
		  
		  for(i=1;i<=brep;i++){
		      idx=mm_sample(msub,.,idinfo,.,1)
			  datastar=dataref[idx,.]
			  Dboot[.,i]=ddfwindow(data0,datastar,k1,k2,gmat,rstype,maxiter,tol)
		  }
		  
		  
		 for(j=1;j<=length(t0uniq);j++){
			  pos=select(1::rows(data0),t0:==t0uniq[j])
			  Kr=length(select(dataref[.,1],t1:<=t0uniq[j]))
			  cires=biasandciddf(D0[pos,1],Dboot[pos,.], k1, cols(data)-k1-2, length(pos), Kr, level,0,msub*sum(t1:<=t0uniq[j]+band & t1:>=t0uniq[j]-band))
			  Dbc[pos,1]=cires.tebc
			  LL[pos,1] =cires.LL
			  UU[pos,1] =cires.UU
			  bias[pos,1]=cires.bias
			  vari[pos,1]=cires.vari
		  }
		  


}





function biasandciddf(  real colvector te,real matrix teboot,  ///
                        real scalar M, real scalar N, real scalar K, real scalar Kr, ///
                        real scalar level,real scalar smoothed, real scalar msub)
{
    if(smoothed){
	    con1=1
		con2=1
	}
	else{
	    con1=msub^(2/(M+N+1))
		con2=Kr^(-2/(M+N+1))
	}
	con3=con1*con2
    bias=J(K,1,.)
	vari=J(K,1,.)
	BovV=J(K,1,.)
	tebc=J(K,1,.)
	UU=J(K,1,.)
	LL=J(K,1,.)
	
	for(i=1;i<=K;i++){
	    TEi=teboot[i,.]
		TEi=select(TEi,TEi:!=.)
		if(length(TEi)<100){
		    bias[i]=.
			vari[i]=.
			BovV[i]=.
			tebc[i]=.
			LL[i]=.
			UU[i]=.
		}
		else{
			
		    bias[i]=con3*(mean(TEi')-te[i])
			vari[i]=variance(TEi')
			BovV[i]=(bias[i])^2 / vari[i] * 3
            tebc[i] = te[i] - bias[i]
            teOverTEhat = (TEi':- te[i] ) * con1	
		    quans =mm_quantile(teOverTEhat,1, ((0.5-level/200)\(0.5+level/200) ))
			//quans* con2
		    LL[i] = te[i] - quans[2] * con2
		    UU[i] = te[i] - quans[1] * con2			
			//quans =mm_quantile(TEi',1, ((0.5-level/200)\(0.5+level/200) ))
			//LL[i] =  quans[1] 
		    //UU[i] =  quans[2] 	
		}
	}
		struct bootdea scalar res
				   res.tebc      = tebc
				   res.vari      = vari
				   res.BovV    = BovV			   
				   res.bias    = bias
				   res.LL 	   = LL
				   res.UU	   = UU
 
        return(res)
}









real scalar function _ddf2( real colvector    X, ///
                            real colvector    Y, ///
						    real colvector    B, ///
						    real matrix       Xref, ///
						    real matrix       Yref, ///
						    real matrix       Bref, ///
						    real colvector    g, ///
					        real scalar       rstype, ///					
					        real scalar       maxiter, ///
					        real scalar       tol)
	{


          k1=length(X)
		  k2=length(Y)
		  M=length(B)+k1+k2
          gX=g[1..k1,1]
          gY=g[k1+1..(k1+k2),1]
          gB=g[(k1+k2+1)..M,1]         
          N=cols(Xref)

	    class LinearProgram scalar q
	    q = LinearProgram()
		c=(1,J(1,N,0))
		lowerbd=.,J(1,length(c)-1,0)
		upperbd=J(1,length(c),.)		
		q.setCoefficients(c)
		q.setBounds(lowerbd, upperbd)

		if(maxiter!=-1){
		  q.setMaxiter(maxiter)
		}
		if (tol!=-1){
		  q.setTol(tol)
		}			

         
  
		Aie1=-gX,Xref
		bie1=X
		Aie2=gY,-Yref
		bie2=-Y
		Aec=-gB,Bref
		bec=B
		if(rstype==0){
			q.setEquality(Aec,bec)		 	
		}
		else{
			lsum=0,J(1,N,1)
			q.setEquality(Aec \ lsum, bec \ 1)		 	
		}
		q.setInequality((Aie1 \ Aie2 ), (bie1 \ bie2))
        beta=q.optimize()		  
         
        return(beta)
		

}




void function _nddf( string scalar     d, ///
                    real scalar       nx, ///
                    real scalar       ny, ///
					string scalar     touse1,  ///
					string scalar     touse2, ///
					string scalar     gname, ///
					string scalar     wname, ///
					string scalar     bname, ///
					real scalar       rstype, ///					
					real scalar       maxiter, ///
					real scalar       tol)
	{

          wmat=st_matrix(wname)
          data=st_data(.,d,touse1)
          data=data'
		  M=rows(data)
          g=st_data(.,gname,touse1)
          g=g'       
          dataref=st_data(.,d,touse2)
          dataref=dataref'
          Xref=dataref[1..nx,.]
          Yref=dataref[nx+1..(nx+ny),.]
          Bref=dataref[(nx+ny+1)..M,.]
          X=data[1..nx,.]
          Y=data[nx+1..(nx+ny),.]
          B=data[(nx+ny+1)..M,.]
 
          theta=J(cols(data),1,.)
  
          for(j=1;j<=cols(data);j++){

                 theta[j]=_nddf2(X[.,j],Y[.,j],B[.,j],Xref,Yref,Bref,g[.,j],wmat,rstype,maxiter,tol)	  
               
         }

          st_view(gen=.,.,bname,touse1)
          gen[.,.]=theta
		

}



real function _nddf2( real colvector     X, ///
                      real colvector     Y, ///
                      real colvector     B, ///
					  real matrix        Xref,  ///
					  real matrix        Yref,  ///
					  real matrix        Bref,  ///
					  real colvector      g, ////
					  real rowvector     wmat, ///
					  real scalar       rstype, ///					
					  real scalar       maxiter, ///
					  real scalar       tol)
	{

          
          nx=length(X)
		  ny=length(Y)
		  nb=length(B)
		  M=nb+nx+ny
          gX=g[1..nx,1]
          gY=g[nx+1..(nx+ny),1]
          gB=g[(nx+ny+1)..M,1]        
          N=cols(Xref)

	    class LinearProgram scalar q
	    q = LinearProgram()
		c=(wmat,J(1,N,0))
		lowerbd=J(1,length(c),0)
		upperbd=J(1,length(c),.)		
		q.setCoefficients(c)
		q.setBounds(lowerbd, upperbd)

		if(maxiter!=-1){
		  q.setMaxiter(maxiter)
		}
		if (tol!=-1){
		  q.setTol(tol)
		}			

          
  

				 Aie1=diag(-gX),J(nx,ny+nb,0),Xref
				 bie1=X
				 Aie2=J(ny,nx,0),diag(gY),J(ny,nb,0),-Yref
				 bie2=-Y
				 Aec=J(nb,nx+ny,0),diag(-gB),Bref
				 bec=B
				 lsum=J(1,length(wmat),0),J(1,N,1)

				if(rstype==0){
				  q.setEquality(Aec,bec)		 	
				}
				else{
					 lsum=J(1,M,0),J(1,N,1)
					 q.setEquality(Aec \ lsum, bec \ 1)		 	
				 }

				 q.setInequality((Aie1 \ Aie2 ), (bie1 \ bie2))

                 beta=q.optimize()		  
                 //q.getInequality()
                 return(beta)
		

}


void function _tenddf( string scalar     d, ///
                    real scalar       nx, ///
                    real scalar       ny, ///
					string scalar     touse1,  ///
					string scalar     touse2, ///
					string scalar     gname, ///
					string scalar     wname, ///
					string scalar     bname, ///
					real scalar       rstype, ///					
					real scalar       maxiter, ///
					real scalar       tol)
	{

          wmat=st_matrix(wname)
          data=st_data(.,d,touse1)
          data=data'
		  M=rows(data)
          g=st_data(.,gname,touse1)
          g=g'       
          dataref=st_data(.,d,touse2)
          dataref=dataref'
          Xref=dataref[1..nx,.]
          Yref=dataref[nx+1..(nx+ny),.]
          Bref=dataref[(nx+ny+1)..M,.]
          X=data[1..nx,.]
          Y=data[nx+1..(nx+ny),.]
          B=data[(nx+ny+1)..M,.]
 
          theta=J(cols(data),rows(data)+1,.)
  
          for(j=1;j<=cols(data);j++){

                 theta[j,.]=_tenddf2(X[.,j],Y[.,j],B[.,j],Xref,Yref,Bref,g[.,j],wmat,rstype,maxiter,tol)	  
                 //q.getInequality()
         }

          st_view(gen=.,.,bname,touse1)
          gen[.,.]=theta
		

}



real function _tenddf2( real colvector     X, ///
                      real colvector     Y, ///
                      real colvector     B, ///
					  real matrix        Xref,  ///
					  real matrix        Yref,  ///
					  real matrix        Bref,  ///
					  real colvector      g, ////
					  real rowvector     wmat, ///
					  real scalar       rstype, ///					
					  real scalar       maxiter, ///
					  real scalar       tol)
	{

          
          nx=length(X)
		  ny=length(Y)
		  nb=length(B)
		  M=nb+nx+ny
          gX=g[1..nx,1]
          gY=g[nx+1..(nx+ny),1]
          gB=g[(nx+ny+1)..M,1]        
          N=cols(Xref)

	    class LinearProgram scalar q
	    q = LinearProgram()
		c=(wmat,J(1,N,0))
		lowerbd=J(1,length(c),0)
		upperbd=J(1,length(c),.)		
		q.setCoefficients(c)
		q.setBounds(lowerbd, upperbd)

		if(maxiter!=-1){
		  q.setMaxiter(maxiter)
		}
		if (tol!=-1){
		  q.setTol(tol)
		}			

          
  

	Aie1=diag(-gX),J(nx,ny+nb,0),Xref
	bie1=X
	Aie2=J(ny,nx,0),diag(gY),J(ny,nb,0),-Yref
	bie2=-Y
	Aec=J(nb,nx+ny,0),diag(-gB),Bref
	bec=B
	lsum=J(1,length(wmat),0),J(1,N,1)

	if(rstype==0){
		 q.setEquality(Aec,bec)		 	
	}
	else{
		 lsum=J(1,M,0),J(1,N,1)
		 q.setEquality(Aec \ lsum, bec \ 1)		 	
	 }

	 q.setInequality((Aie1 \ Aie2 ), (bie1 \ bie2))

	 dval=q.optimize()

	 if (q.converged()==1){
	 	bslk=q.parameters()
	 	beta=dval,bslk[1..M]
	 }
	 else{
	 	beta=J(1,1+M,.)
	 }
              	  
     //q.getInequality()
     return(beta)
		

}






mata mlib create lgtfpch, replace 
mata mlib add lgtfpch *()
mata mlib index

end